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ABSTRACT 

Stoichiometric constraints play a role in the dynamics of natural populations, but 
are not explicitly considered in most mathematical models. Recent theoretical works 
suggest that these constraints can have a significant impact and should not be neglected. 
However, it is not yet resolved how stoichiometry should be integrated in population 
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dynamical models, as different modeling approaches are found to yield qualitatively dif- 
ferent results. Here we investigate a unifying framework that reveals the differences and 
commonalities between previously proposed models for producer-grazer systems. Our 
analysis reveals that stoichiometric constraints affect the dynamics mainly by increasing 
the intraspecific competition between producers and by introducing a variable biomass 
conversion efficiency. The intraspecific competition has a strongly stabilizing effect on 
the system, whereas the variable conversion efficiency resulting from a variable food 
quality is the main determinant for the nature of the instability once destabilization 
occurs. Only if the food quality is high an oscillatory instability, as in the classical 
paradox of enrichment, can occur. While the generalized model reveals that the generic 
insights remain valid in a large class of models, we show that other details such as 
the specific sequence of bifurcations encountered in enrichment scenarios can depend 
sensitively on assumptions made in modeling stoichiometric constraints. 
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1. Introduction 



Many ecological models quantify energy and biomass flow solely in terms of carbon, whereas 
stoichiometric constraints, arising in part from different nutrient ratios in the populations, are 
only captured indirectly. Recently, it has been shown that already minor extensions can make 
carbon-based models stoichiometrically explicit and thus sign ificantly enhance the qu alitative un- 
derstanding of laboratory experiments and field observations (jSterner and Elseril2002l ). 



Stoichiometric constraints have been found to affect particularly the conversion efficiency 
from the first to the second trophic level of an ecosystem and the rate of primary production 
(| Andersen et al.ll2004l ). To understand the reason for the variable conversion efficiency, consider 
that most primary producers are flexible in their use of nutrients and are thus characterized by 
highly variable nutrient content. By contrast, grazers have a relatively fixed internal stoichiometry. 
Thus, not all carbon available to the grazer can be utilized if the nutrient concentration in the 
producer biomass is low. 

The produce r's nutrient content depends on many complex processes governing nutrient flows 
(jPeAngelisI 1 19921 ). but is particular ly depen dent on grazing. Although grazing can enhance the 
recycling of nutrients in the system (ISternerill98& ). it also sustains accumulation of nutrients in the 
biomass of the higher trophic levels. In particular in systems in which the recycling of nutrients is 
essential, the accumulation of nutrients in the grazers can lead to a depletion of available nutrients. 



- 3 - 



Thereby biomass in higher trophic levels can affect primary production, the nutrient content of the 
producers and consequently also the conversion efficiency of the grazer. 

Nutrient accumulation and variable conversion efficiency introduce a complex feedback mech- 
anism because the rate of primary production and the growth of grazers become dependent on the 
biomasses of all populations in th e system. Thus, stoichiornetric mechanisms can allow fo r coexis- 



tence of competing grazer species (|Hallll2004l : lLoladze et al.ll2004l : iDeng and Loladzd 120071 ) and can 



prevent local extinction of producer species in heterogeneous habitats (jMiller et al.ll2004l ). Even in 



affect the system dynamics ( 


Huxel 


1999; 


Loladze et al. 


2000: 


Muller et al. 


2001; 


Sterner and Elser 


2002; 


Grover 


2003: 


Kooiiman et al. 


2004 


)• 



A point of particular concern is that seemingly similar models can exhibit different dynamics, 
depending on the functional forms that are used to describe the conversion efficiency. For instance, 
in some models stable stationary behavior even at high levels of enrichment and multi-stability are 
possible whereas in others enrichment leads to deterministic extinction. Since the metabolism of 
even a single cell is highly complex, every specific mathematical function, formulated to describe 
stoichiometric constraints on the level of the population, necessarily involves strong assumptions. It 
is therefore an important practical challenge to identify the decisive feature of the functional forms 
that determine the dynamics of the populations and the producer-grazer system and therefore have 
to be captured in order to formulate credible ecosystem models. 



^In this paper we use the approach of generalized modeling ([Gross and Feudelll2006l ; iGross et al 

20091 ) to analyze the effects of stoichiometric constraints. In a generalized model the rates of 



processes do not need to be restricted to specific functional forms, allowing for the investigation of 
a large class of models. We then compare the results of the generalized model to several specific 
examples. 

We show that the generalized model provides a unifying framework that explains differences 
and commonalities between the different specific models, whereas the specific models allow for 
a detailed numerical investigation revealing additional insights. Our analysis of the generalized 
stoichiometric model identifies six parameters that capture the stability of all stationary solutions 
in a large class of different models. By combining the generalized analysis with the investigation of 
specific models we then show that intra-specific competition in the producer is the main determinant 
of stability. When stability is lost, the decisive factor determining the nature of the instability is 
the variability of the biomass conversion efficiency. If the conversion efficiency is approximately 
constant then the destabilization is caused by an oscillatory instability as in the classical paradox 
of enrichment (|Rosenzweig)ll97ll ). If however the conversion efficiency becomes strongly dependent 
on nutrient content of the producer then the oscillatory instability is replaced by a different type 
of instability that can be related to the onset of an Allee effect. Our generalized analysis thereby 
provides a deeper understanding of the dynamics observed in specific models and shows that the 
avoidance of the paradox of enrichment and the appearance of Allee effects are generic features of 
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stoichiometric models. 



2. A generalized food chain model with variable efficiency 

We consider the class of predator-prey models in which a primary producer X is consumed by a 
grazer Y. To conform with the specific models studied in Sec. [4] let X and Y measure the respective 
biomass densities in terms of carbon concentrations. We assume that grazing is the only source of 
mortality for the producer, whereas the consumer suffers from natural mortality at a constant rate 
D. The dynamics of X and Y can then be expressed as 



dt' 



-X 



S{X,Y)- F{X)Y 



—Y = E(X, Y)F(X)Y - DY . 
at 



(1) 
(2) 



where S{X,Y), E{X,Y), and F{X) denote the primary production, the conversion efficiency and 
the functional response of the grazing, respectively. 

Typically, the first step in the analysis of a model, such as Eqs. ([l]l2]), is to restrict the func- 
tions to specific functional forms, for instance by assumin g that the production follows logistic 



growth. An alternative route, sometimes taken in ecologv (iGardner arid Ashbvlll97Cl: iMav 



Murdoch and OatenI Il975l : iDeAngelis 



1975 



LevinI 119771 : Murdoch 



197' 



Wollkind et al 



1972 



[1982), is 



to analyze the whole class of models described by the set of equations containing the general, 
i.e., unrestricted, functi ons. Here we analyze the general produce r- grazer model by the approach 



of generalized modeling ()Gross and Feudel 



2006 



Gross et al.ll2009l ). which focuses on the dynamics 



close to positive, i.e., feasible, steady states. In such a state the right-hand side of Eqs. ([HIS]) van- 
ishes, so that the system can reside in the steady state infinitely long. A steady state can be stable, 
corresponding to an equilibrium of the system, or unstable, so that even a small perturbation leads 
to departure form the vicinity of the steady state. The aim of generalized modeling procedure is to 
capture the decisive factors governing the stability of steady states in the unrestricted generalized 
model by a set of parameters having clear interpretations. We present an outline of this procedure 
below, whereas a detailed description is provided in the supporting material lA.ll 

The biomasses of the consumer and grazer populations in the steady state generally depend on 
many details of the functional forms in the model and hence cannot be computed in the generalized 
model. However, we can formally denote the biomasses in an unknown steady state as X* and 
Y*, respectively. This allows us to map the unknown steady state to a defined value by a suitable 
normalization. Below, we use lower case letters to denote variables {x,y) and functions (e, /, s) that 
have been normalized to unity in the steady state under consideration, e.g. f{x) = F{X) / F{X*). 
We emphasize, that for feasible steady states this normalization is possible without further as- 
sumptions. Specifically, no further assumptions are made on the number of steady states, or the 
value of biomasses in the normalized steady state. Further, we set the biomass turnover-rate of the 
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producer to one by means of a timescale normalization, which is equivalent to measuring all rates 
in multiples of the producer's biomass turnover. We thus obtain 



dt 
d_ 

dt 



X 



y 



six,y) - f{x)y , 
r{e{x,y)f{x)y - y) 



(3) 
(4) 



where the new quantity r is defined as the turnover-rate of the grazer in the steady state. By the 
normalization we have thus moved from a model that did not contain any parameters to a model 
containing the unknown parameter r. As the parameter is measured in multiples of the produce r 
turnover-rate it can be assumed that its value is between and 1 for most systems (|Hendrikall999l ). 



The stabilit y of steady states to small pert urbations can be computed from the so-called 
Jacobian matrix (jGuckenheimer and Holmesll2002l ). which in the present model can be written as 



dx - 7 cr^ 
riVx + 1) 



y 

rr]y 



1 



(5) 



where 



7 


= imi 






ay 






= m<^^y) 







x=x*,y=y* 

x=x*,y=y* 
x=x* ,y=y* 

x=x* ,y=y* 



(6) 



Equation ([6]) defines constant scalar quantities that can be treated as unknown parameters. To 
understand the meaning of these parameters consider their relation to the original function, e.g. 



O'x 



§^s{x,y) 



x=x* ,y=y* 



(7) 



X=X-'.Y=Y* 



If S were any linear function, then the parameter ax were equal. More generally every power- 
law, such as S{X) ~ AXP, corresponds to ax = p- Even if S is a more general function the 
corresponding p arameter ax deiiotes t he so-called elasticity, a nonlinear measure for the sensitivity 
of the function (jFell and Saurd Il985l ) . Let us emphasize that the elasticities ca n be estimated 
direc tly from field observations by means of conventional linear regressions (e.g., (jQian and Giles 
20071 )). In contrast to parameters defined for mathematical convenience, such as half-saturation 
constants, the generalized parameters do not require reference to an artificial state (e.g., the half 
saturation point), which may be far from the natural state of the system. 

In a system of differential equations the stability o f steady states de pends on the eigenvalues 
of the Jacobian, which are in general complex numbers (jKuznetsovl |2004| ) . A steady state is stable 
whenever all eigenvalues of the Jacobian have negative real parts. When parameters are changed 
the stability of a steady state is lost if the change causes at least one eigenvalue to acquire a positive 
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real part. If a steady state becomes unstable small fluctuations in the population densities cause 
the system to depart from this equilibrium. The dynamical transition involved in a loss of stability 
is called bifurcation whereas the set of critical parameter values at which the transition occurs is 
called the bifurcation point. 

For real matrices, such as J, the loss of stability can occur in two generic scenarios: a) In 
a Hopf bifurcation the stability of the steady state is lost and either a stable limit cycle emerges 
(supercritical Hopf) or an unstable limit cycle vanishes (subcritical Hopf). A Hopf bifurcation 
therefore marks a transition between stationary and oscillatory population densities, where the 
oscillations are transient in the subcritical case and sustained in the supercritical case, b) In a 
saddle-node bifurcation the steady state under consideration collides with another steady state. In 
general the steady states annihilate each other, so that the system approaches some other attractor. 
This transition is in general not reversible because the system will typically not return to the original 
state even if its stability is restored by subsequent changes of parameters. The bifurcation therefore 
often marks the onset of an Allee effect (e.g. see Seed]). Because of certain symmetries a degenerate 
form of the saddle-node bifurcation, the transcritical bifurcation is sometimes encountered. In this 
bifurcation the two steady states cross, exchanging their stability. In ecological models this soft 
transition is often encountered when the grazer population in a steady state becomes positive. 

In the present model a saddle-node bifurcation occurs when 

Vyi(^x - 7) - {(^y - l)(??x + 7) = (8) 

and a Hopf bifurcation when 

(Tx - 7 + rr]y = (9) 

and 

Vy{(^x - 7) - {cTy - l)(??x + 7) > 0. (10) 
The computation of these conditions is shown in the supporting material IA.2[ 

By means of the normalization procedure we have managed to identify the decisive properties 
that determine the stability of all steady states in all models of the form of Eqs. ([T]l2]). In Sec. [3] we 
map these parameters (r, 7, ax, o'y, rjx, rjy) to a rescaled parameter set (r, 7, c^, Cy, n^, Uy) that 
allows for a clear ecological discussion. 



3. Generalized analysis 



For exploring the dynamics of the generalized stoichiometric model our main tool will be 
three-parameter bifurcation diagrams, such as the one shown in Fig. 1. These d i agrani s can be 
constructed from the bifurcation conditions Eqs. (fHl fTO]) by following IStiefs et al.l (j2008l ). Every 
point in the diagram corresponds to a steady state that is characterized by a specific combination 
of the three parameters on the axis. In the three-dimensional space the bifurcation points form 
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surfaces, separating volumes with different local dynamics. Except for diagrams shown in top- 
view, we have oriented the diagrams such that the steady states in the topmost volume are stable. 
Changing the parameters causes a loss of stability if and only if it leads to a departure from the 
topmost volume by crossing one of the surfaces. Throughout we use red surfaces to indicate Hopf 
bifurcations, marking the onset of oscillations and blue surfaces to indicate saddle-node bifurcations, 
which can be related to the onset of an Allee effect (see Sec. [4]). Additionally, there is one more 
surface, corresponding to transcritical bifurcations, which is not shown because it is located in the 
front plane of the diagram. On this surface the biomass of the grazer vanishes in the steady state, 
corresponding to an extinction of the grazer population. 

From previous studies it is known that the relative turnover-rate of the consumer has relatively 
little impact on the dynamics. This is also confirmed by our equations for the bifurcation surfaces. 
The location of the saddle-node bifurcation surface is independent of r and the Hopf bifurcation 
surface is affected only mildly. In the following we set this parameter to a moderate value of 0.3. 

A very important parameter for the dynamical stability is 7, the sensitivity of the grazing 
rate on the biomass of the producer. This parameter is 1 if the grazing rate depends linearly on 
the producer, e.g., in a Lotka-Volterra model. The parameter can be as low as 0, if the grazing 
saturates, or as high as 2 if the grazing rate depends quadratically on the producer. The latter 
case is encountered for instance in models with a type-III functional response in the limit of scarce 
prey. 



Gross et al.l (|2004l ) showed that the parameter 7 has a strong stabilizing influence on the 
system. Increasing the nutrient or energy supply to the system can therefore lead to a destabilization 
of the system, because the increasing number of producers leads to increased saturation of the 
grazer and hence lowers 7. This mechanism is closely linked to Rosenzweig's paradox of enrichment 



(jRosenzweig and MacArthur 



1963 



Rosenzweig) Il97ll ). which specifically describes destabilization 



in a Hopf bifurcation. Originally, it was felt paradoxical that nutrient or energy enrichment, while 
supposedly beneficial for the individuals, could lead to destabilization of steady states and an 
increased risk of extinction. From a modern perspective the paradoxical aspect of Rosenzweig's 
mechanism is rather that it is found quite generally in models, but is only rarely observed in nature 
(jMorin and Lawlerill995l ). 



Among several other explanations for the absence of the paradox of enrichment in nature, it has 
been noted that the destabilizat ion in a Hopf bifurcation can be avoided in stoichiometric models 



( Loladze et al 



2000 



Diehll 120071 ). Specifically, the absence of the Hopf bifurcation can be linked to 



the variable conversion efficiency, captured in our model by the parameters r]x and rjy. If nutrient 
content of the producer is not limiting then the conversion efficiency is constant and both parameters 
are zero. If, due to nutrient limitation, the conversion efficiency decreases with increasing densities 
of the producer or consumer, the respective parameter r]x or r]y becomes negative. In general both 
parameters are closely linked because both scale inversely with the value of conversion efficiency. In 
the extreme limit when the variable conversion efficiency becomes zero both parameters approach 
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Fig. 1. — Bifurcation diagrams of a generalized producer-grazer model. At constant conversion 
efficiency, = Uy = 1, a Hopf bifurcation (red) is the only source of instability. In models with 
variable conversion efficiency the Hopf bifurcation is replaced by a saddle-node bifurcation (blue) 
if food quality is low, leading to an avoidance of the paradox of enrichment. Every point in the 
diagram corresponds to a steady state in specific models. These states are stable if the point is 
located above all bifurcation surfaces and unstable otherwise. Parameters: r = 0.3, c„ = 0. 
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minus infinity. It is therefore advantageous to introduce rescaled parameters = 1/(1 — r]x) and 
Uy = 1/(1 — r]y) where 1 now corresponds to a constant conversion efficiency and to a vanishing 
conversion efficiency. These parameters can be interpreted as direct indicators of the grazer's food 
quahty, ranging from optimal (1) to worthless (0). 

The effect of food quality on the dynamics of the general stoichiometric model is shown in Fig. 1. 
In the figure we have set Ux and Uy to identical values, which can be justified in specific models 
and does not affect our conclusions (see Sec. [Hand supporting material IA.3p . If the conversion 
efficiency is constant (tLx = Uy = 1) then saddle-node bifurcations cannot occur for positive values 
of the grazer's sensitivity 7. Hence, the Hopf bifurcation is the only remaining source of instability 
in feasible models (see supporting material I A. 41 for a mathematical proof). However, if the food 
quality is decreased, a saddle-node bifurcation can occur. The surface of Hopf bifurcations ends 
when it meets the saddle-node bifurcation surface. The saddle-node bifurcation thus replaces the 
Hopf bifurcation as the primary source of instability. In particular Hopf bifurcations cannot occur 
if the food quality is below 0.5. 

Our results on the effect of variable conversion efficiency confirm and generalize the previous 
observation that the paradox of enrichment, i.e., the onset of oscillations, can be avoided in stoi- 
chiometric models. However, it also clearly shows that the instability is not avoided entirely, but 
rather that oscillations are replaced by a different type of instability that may lead more directly 
to extinction. 

Let us now focus on the final two parameters ax and ay, which denote the sensitivity of produc- 
tion to the producer and grazer populations, respectively. If producers are not limited by nutrients 
or other factors then it is reasonable to assume that production is linear in the number of producers 
and independent of the number of grazers {ax = 1 and ay = 0). Lower values of the parameters are 
found if production suffers from the sequestration of nutrients in the producer or grazer population. 
Again, in the extreme limit where sequestration precludes further production both of the parame- 
ters go to minus infinity. We therefore discuss the impact of nutrient sequestration in terms of the 
rescaled parameters Cx = (1 — ax)/ (2 — ax) and Cy = —ay /{I — ay). The parameters Cx and Cy can 
be interpreted as measures of the intraspecific competition for nutrients in the producer population 
and the interspecific competition for nutrients between producer and grazer, respectively. If Cx is 
{cy is 0) then the production by a given individual is independent of the density of other producers 
(grazers). By contrast, a value of 1 indicates that all accessible nutrients have been taken up so 
that no further production is possible. 

For the stability, the interspecific competition, Cy, is of lesser importance. From the bifurcation 
conditions, Eq. ([9]), it is apparent that the Hopf bifurcation does not depend on this parameter. 
As we have seen that the Hopf bifurcation is the only source of instability in models with constant 
conversion efficiency the interspecific competition cannot affect the stability in such models. In 
models with variable conversion efficiency the parameter affects the saddle-node bifurcation, but 
its impact is relatively mild. 
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It can be seen from Fig. 1, increasing the intraspecific competition for nutrients, Cx, has a 
strongly stabihzing effect. As the parameter is increased the bifurcation surfaces drop away, so 
that the stable topmost volume is enlarged. In this sense the effect of changes in the intraspecific 
competition is reminiscent of the paradox of enrichment. We therefore use the term paradox of 
competition to denote the observation that strong intraspecific competition lends stability, while 
decreasing intraspecific competition harbors instability. Note however that this paradox is partially 
resolved if one takes into account that increasing competition also brings the system closer to the 
= 1 plane, where the primary production vanishes and the grazer goes extinct in a transcritical 
bifurcation (see supporting material IA.5p . Although no loss of stability is involved in this case, it 
shows that increasing competition can after all have a negative impact on the system. 

We summarize the results of the general analysis in three basic statements: First, variable 
conversion efficiency can have an important impact on the dynamics. In particular the oscillatory 
instability (Hopf bifurcation) is replaced by a different type of instability (saddle-node bifurcation) 
when variable conversion efficiency is introduced. This leads to an avoidance of the classical para- 
dox of enrichment, but does not convey increased stability. Second, increasing the intraspecific 
competition has a stabilizing effect comparable to increasing the sensitivity of the grazer, giving 
rise to a paradox of competition. Third, the interspecific competition for nutrients, i.e., the seques- 
tration of nutrients in higher trophic levels, has little effect on the stability apart from potentially 
increasing the intraspecific competition for the remaining accessible nutrients. 



4. Specific stoichiometric modeling approaches 



A more detailed understanding of the dynamics can be gained if the generalized analysis is com- 



bined with t he investigation of specific models. We consider a mode l proposed bvlKooijman et al.. 
(2004, 2007 ) based on the dynamic energy budget (DE B) theory (|Kooiimanl I2OI0I ). a simplifiec. 



DEB model (SDEB), the model by (jLoladze et al.ll2000l ) (LKE), which is based on Liebig ' s min - 
imum law, a smooth analogon (SA) to the LKE model, and the related model by iDiehll (j2007l ) 
which considers the effect of light limitation (DLL). We emphasize that these models differ by the 
enrichment scenario they consider as well as the mathematical representation. 

The DEB model captures that variable nutrient content of producers mainly arises from storage 
compartments for nutrients. The nutrients in the producer's structure and in storage are therefore 
represented separately in these models. It is assumed that the total amount of nutrients in the 
system is fixed, and that there are essentially no free nutrients because of the fast uptake kinetics of 
the producers. Under these assumptions the conservation law for nutrients can be used to express 
the nutrients in the producer's storage as a function of the nutrients in the producer's structure and 
the nutrients in the grazer, leaving the latter two as the only dynamical variables in the system. 
Therefore, a system of the form of Eqs. ([T]l2|) is obtained. The specific rate laws in this system can 
be derived using the concept of synthesizing units, which represents th e population as individual 



entities assembling biomass from distinct substrates (iMuller et al 



2001 



)• 
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Although the nutrient-to-carbon ratio in the producer's structure is fixed, the food quahty 
in the DEB model is variable because the producers carry a variable amount of nutrients in their 
storage. A significant simplification is obtained in the SDEB model where it is assumed that the 
grazer cannot utilize nutrients from the producer's storage compartment, fixing the conversion 
efficiency. 

The LKE model, while very similar to the DEB model, represents the nutrients in the producer 
differently. Instead of distinguishing between structure and storage, the model explicitly considers 
the amount of carbon and phosphorus in the producer. Assuming a fixed amount of phosphorus in 
the system and a fixed carbon-to-phosphorus ratio in the grazer, the total amount of phosphorous 
that is accessible to the producer is given by a conservation law for phosphorus. Thus a model 
of the form of Eqs. ([l]l2]) is obtained where now the carbon in the producer and the grazer are 
the dynamical variables, while the phosporus in the prey is given by the conservation law. The 
expressions governing nutrient uptake and coversion efficiency are based on Liebig's minimum law 
and therefore depend either on carbon or phosporus depending on which is most strongly limiting. 
This implies that the conversion efficiency is constant when the grazer is carbon-limited but depends 
on the phosporus content of the producer if the grazer is phosphorus- limited. 

The use of Liebeigs law, which induces non-smooth behavior in the LKE model, is well justified 
when one factor imposes a much stronger limitation than others in the system. However, if limiting 
factors are of comparable importance microscopic fluctuations on the level of the individuals may 
lead to transient shortages in the second most limiting factor. In the SA model we captured this 
effect by again applying the concept of synthesizing units. We thus obtain a model in which the 
conversion efficiency is approximately constant when carbon is limiting, but smoothly decreases 
with increasing importance of phosphorus. 

The closely related DLL model describes an aquatic system that is limited by nutrients and 
energy, where the energy input through light is controlled by a parameter that describes the depth 
of the mixed layer. It is assumed that the two constraints act multiplicatively, leading to a smooth 
transition between light-limitation and nutrient-limitation. The grazer growth limitation by carbon 
and nutrients is modeled by the synthesizing units approach. The full model contains, besides the 
producer and grazer b iomasses, two additional variables describing dissolved and detrital nutrients. 



However, iDiehll (j2007l ) shows that the dynamics of the full model can already be observed in a 
simplified model where the detrital nutrients are set to zero and the dissolved nutrients are fixed by 
a conservation relation. In the model variable conversion efficiency enters through a factor Q which 
describes the nutrient content per producer carbon biomass. When the availability of nutrients in 
the system increases Q, increases up to a maximal value at which the producers stop assimilating 
further nutrients. The conversion efficiency is therefore constant if the nutrient availability is high, 
but becomes variable once the nutrient content of producers falls below the maximal value. 

The dynamics of the specific models are depict ed in the 1-parameter diagrams shown in Figs. 2 - 



4, which we reproduced using the software AUTO (jDoedel et al.lll997l : lDoedel and Oldemanll2009l ). 
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These diagrams show the grazer biomass in the steady state, as a function of one of the models 
parameters. Sohd hnes in the diagrams correspond to stable steady states, whereas dashed lines 
correspond to unstable steady states. If limit cycles exist, the upper and lower turning points of 
the cycle are indicated by thin lines, which are solid if the cycle is stable and dotted if the cycle is 
unstable. In all diagrams we have chosen the parameter being varied such that, if read from left 
to right, the diagram represents an enrichment scenario. The respective parameters are the total 
amount of nutrients in the SDEB and DEB models, the producer carrying capacity in the LKE 
model, the total amount of carbon in the SA model, and the depth of the mixed layer in the DLL 
model. 

Although the models are closely related some important differences exist between the bifur- 
cation diagrams. In all models there is a threshold below which the grazer population cannot be 
sustained. If this threshold is exceeded the system is sufficiently enriched to support grazers. In 
four of the five models (SDEB, LKE, SA, DLL) this lower threshold corresponds to a transcritical 
bifurcation (TC), representing a soft transition. If the enrichment parameter is decreased toward 
the threshold then the grazer population declines gradually until it vanishes in the bifurcation 
point. When the enrichment parameter is subsequently increased again, the grazer can re-invade 
the system immediately beyond the bifurcation point because the steady state with zero grazer 
density becomes unstable in the transcritical bifurcation. By contrast, in the DEB model, the 
enrichment threshold corresponds to a saddle-node bifurcation representing a hard transition. If 
the enrichment parameter is lowered below the threshold in this model the positive steady state is 
annihilated causing a rapid decline of the grazer population. Even if the enrichment parameter is 
subsequently increased beyond the threshold, the grazer cannot re-invade the system immediately 



as the steady state at zero grazer density remains stable (IScheffer et al.ll200ll ). 



If the enrichment parameter is sufficiently increased then all models undergo a supercritical 
Hopf (H) bifurcation giving rise to predator-prey oscillations. However, only in Fig. 2 (top right) 
corresponding to the SDEB model these oscillations persist even for very high levels of enrichment. 
In the four other models the oscillations end when the stable limit cycle is destroyed in subsequent 
bifurcations. In the DEB, LKE, and SA models the stable limit cycle vanishes as it collides with 
an unstable steady state in a global homoclinic bifurcation. Although a detailed discussion of this 
type of bifurcation exceeds the scope of this paper, let us mention that the homoclinic bifurcation 
can cause repeated outbreaks and, in higher dimensional systems, chaotic dynamics. By contrast, 
in the DLL model the stable limit cycle is destroyed as it collides with an unstable limit cycle in 
a so-called fold bifurcation (F) of cycles. Although this is another example of a hard transition, 
in the DLL model, this transition occurs very close to a stable steady state in which the system 
subsequently settles. Therefore, in practice this transition is hard to distinguish from a smooth 
cessation of oscillations. 

Finally, for very high levels of enrichment an upper threshold is encountered in two of the 
models (LKE, SA). Beyond this transition the nutrient content of the producer is too low to meet 
the demands of the grazer, so that the grazer population vanishes in another soft transition. In 
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the SDEB model the grazer can persist in an oscillatory state. However, for strong enrichment the 
amplitude of the oscillations grows so that resilience is low and extinction becomes likely. In the 
DEB model the extinction occurs at the homoclinic bifurcation in a hard transition. By contrast, 
the DLL model, at least for the present parameter set, can sustain a stationary level of grazers at 
high levels of enrichment. 

For understanding the differences and commonalities between the modeling approaches it is 
conductive to relate the results from the specific model to our generalized analysis. Although all 
of the specific models are parameterized differently the stability of every single steady state in all 
of the models is captured by the generalized analysis and can be expressed as a function of the 
generalized parameters. We could therefore fix the parameters in one of the specific models, pick 
a steady state, and map this particular steady state to a point in the six-dimensional parameter 
space of the generalized model. When a parameter of the specific model is changed smoothly, 
the steady state under consideration will move on a path through the parameter space of the 
generalized model. In the center and lower panels of Figs. 2-4 we show such paths corresponding 
to the enrichment scenarios studied above. In general, all six generalized parameters can vary as 
one of the specific parameters is changed. We have therefore varied the food quality parameter 
Uy according to the specific value that they locally assume in the respective steady state, whereas 
we set the remaining two parameters, the grazer turnover-rate r and interspecific competition Cy 
to fixed values. Comparing the three-parameter bifurcation diagrams reveals that changing these 
parameters does not affect the bifurcation diagram qualitatively. 

The simplest behavior is found in the SDEB model. In the model the conversion efficiency is 
constant so that the parameters describing the food quality in the generalized model remains at 
1. As the nutrient content of the prey is not important the grazer population can be sustained 

whenever production takes place. As noted above the only real source of instability is the Hopf 
bifurcation. Therefore, the bifurcation diagram in this and other constant efficiency models docs not 
contain further complications beyond what is known from classical models in which stoichiometric 
constraints are neglected, e.g., the Rosenzweig-MacArthur model. 

In the DEB model, shown in Fig. 2 (right), the variable conversion efficiency has important 
consequences. If the nutrient level is low then the system can undergo a saddle-node bifurcation 
leading to a strong Allee effect. In this case the limiting factor for the grazer is not the density 
of producers, but rather the total amount of nutrients that can be harvested. Therefore, a mini- 
mum density of grazers is necessary to control the number of producers and thereby maintain the 
nutrient-content per producer at a high level. If the grazer density falls below this level (roughly 
corresponding to the unstable steady state) the producer escapes the grazer's control and the nu- 
trient content per producer becomes so low that even when the grazing rate saturates, not enough 
nutrients can be harvested to sustain the grazer population. This risk is particularly pronounced if 
the producer and grazer are both strongly nutrient limited. Then a decreasing number of grazers 
will lead to an increased number of producers and therefore to a further dilution of the nutrients 
in the producer population. This constitutes the positive feedback causing the destabilization in 
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SDEB DEB 




Fig. 2. — Bifurcations diagrams of the DEB model (right) and the SDEB model (left) for nutrient 
enrichment. In the one-parameter diagram wide lines mark the poistion of steady states, which can 
correspond to stable equilibria (solid) or unstable states (dashed). The this solid lines show the 
upper and lower turning points of a limit cycle. Letters mark the transcritical (TC), saddle-node 
(5) and Hopf (H) bifurcation points of steady states and a homoclinic bifurcation G, in which 
the limit cycle is destroyed. The one parameter bifurcation diagram corresponds to a path in the 
bifurcation diagram of the generalized model, which is shown from two perspectives in the center 
and lower panels. The bifurcation points in the specific models correspond to crossings of Hopf (red) 
and saddle-node (blue) bifurcation surfaces in the generalized model. A transcritical bifurcation 
surface is located in the front plane of the generalized model and is hence not shown. 
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Fig. 3. — Bifurcation diagrams Bifurcations diagram of the LKE (left) and the SA model (right) 
in analogy to Fig. [2j In contrast to the bifurcations observed in the SDEB and DEB model two 
disconnected branches of solutions exist in the bifurcation diagrams for the LKE model, which have 
been indicated by lines of different color. The separate branches appear due to the non-smooth 
behavior of Liebig's minimum that was used in the model. This non-smoothness is removed in the 
SA model. 
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Fig. 4. — Bifurcation diagram of the DLL model (top) in analogy to Fig. [2j Below the additional 
unstable steady state (gray) the grazers go extinct which leads to an Allee effect. The DLL model 
behaves more smoothly than the LKE model, but still contains a minimum function which causes 
a jump (T) in the generalized bifurcation diagrams (center and lower panel). Although the shape 
of the path that the model takes is similar to the SA model, the model solution branch exits the 
unstable region through the Hopf bifurcation surface instead of the saddle-node bifurcation surface, 
causing the appearance of a second Hopf bifurcation H2 in the one-parameter bifurcation diagram 
of the specific model (top panel). At a fold bifurcation {F) the stable limit cycle collides with an 
unstable limit cycle. 
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the saddle- node bifurcation and the appearance of the Ahee effect. 

Another, more subtle consequence of the variable conversion efficiency is the homoclinic bi- 
furcation. Although we cannot compute the surface of homoclinic bifurcations in the generalized 
model, the existence of such a surface is already evident from the fact that the Hopf bifurcation 
surface ends as it meets the saddle-node bifurcation surface. This is explained in more detail in the 
supporting material IA.6[ 

The LKE model, shown in Fig. 3 (left), contains two distinct solution branches. In one of 
the system is limited by carbon so that the food quality remains at one. This branch exhibits the 
same dynamics as the SDEB model. By contrast, the steady state on the other branch the system 
is phosphorus- limited. For high levels of enrichment this branch vanishes in a soft transcritical 
bifurcation, if the phosphorus content of the producer becomes to low to satisfy the grazer's demand. 
While the transition is reminiscent to the saddle-node bifurcation in the DEB model, the two 
transitions take place if enrichment exceeds a threshold in the LKE model or falls below a threshold 
in the DEB model. This difference arises because different enrichment scenarios have been studied. 
In contrast to the DEB model, where the total amount of available nutrients was increased, the 
increasing carrying capacity for the prey in the LKE model can lower the prey's nutrient content. 

The existence of two solution branches in the LKE model can be related to the discontinuous 
nature of the minimum law. This is revealed by the SA model, shown in Fig. 3 (right), where 
the two solution branches connect in the shape of a parabola. However, the topological differences 
between LKE and SA occur in the space below the bifurcation surfaces where they only affect the 
unstable steady states. The stable solution branches, which are of primary importance for ecology, 
are therefore similar in the two models, whereas the structure of the unstable branches. 

The DLL model, shown in Fig. 4 also contains two different regimes corresponding to light 
limitation and nutrient limitation. The DLL model is smoother than the LKE model in the sense 
that the two regimes do not correspond to distinct solution branches in the bifurcation diagram. 
However, in the generalized bifurcation diagram the non-smoothness arising from the nutrient 
saturation of producers is still discernible as a sudden jump (T). At low levels of energy input, the 
nutrient level is above the saturation threshold. Therefore, the conversion efficiency is constant 
and the DLL system behaves as the SDEB model. If the availability of energy increases such 
that nutrient limitation sets in the food quality starts to drop. In contrast to the SA model the 
solution branch does not leave the unstable volume through the saddle-node bifurcation surface 
but exits through the Hopf bifurcation surface. Therefore the return to stationarity at high level of 
enrichment occurs by the Hopf and fold bifurcations rather than the saddle-node and homoclinic 
bifurcations observed in the LKE and SA models. 

The comparison from the previous paragraph suggests that the dynamics in the nutrient- 
limited regime of the DLL model is more closely related to the DEB model than it is to the 
respective branches in the LKE and SA models, although the branch is followed in the opposite 
direction because energy supply rather than nutrients is increased in the DLL model. Comparing 
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the DEB and DLL models one notes the different bifurcations connecting stationarity and oscillatory 
behavior, mentioned above. More importantly, the saddle-node bifurcation, which causes the strong 
Allee effect in the DEB model seems to be absent in the DLL model. However, a minor decrease 
of the total nutrient concentration is sufficient to observe the saddle-node bifurcations in the DLL 
model. In this case also the sequence of bifurcations from the DEB model reappears in the DLL 
model. 



5. Summary and Conclusions 

In the present paper a generalized stoichiometric model was analyzed and compared to five 
specific models. The generalized analysis showed that the stability of steady states depends on six 
parameters with clear ecological interpretations. Among these parameters the analysis identified 
the indicators of food quality and intraspecific competition in the producer as having a particularly 
strong impact on the dynamics. Specifically, we found that intraspecific competition is the main 
determinant of stability in the system while food quality determines the nature of the instability 
once destabilization occurs. 

By connecting observations that have previously been made in isolation in specific models, 
generalized modeling can facilitate the comprehension of the bigger picture that persists irrespective 
of specific decisions made in the modeling process. Instead of repeating the detailed discussions 
from the main part of the paper, let us briefly summarize the generalized picture and discussing 
the generic enrichment scenarios that can be expected in stoichiometric models. We start out 
with the situation in which producers are limited by external factors such as light. In this case 
the intraspecific competition between producers is strong and the nutrient content per producer is 
high. In this case the grazer population is clearly limited by the density of producers. 

As the system is enriched by relaxing the external constraints a threshold at which the density 
of producers is high enough to meet the essential carbon and energy demands of grazers. This 
threshold corresponds in general to a transcritical bifurcation in which grazers enter the system 
smoothly. Once the grazer is established in the system the grazer population will control the density 
of producers. As the system is further enriched the intraspecfic competition between consumers 
is relaxed while the nutrient content of consumers, i.e., the food quality, remains high. In this 
case the destabilization is well captured by models that are not stoichiometrically explicit or use 
constant biomass conversion efficiencies. In these models a supercritical Hopf bifurcation causing 
a transition to predator-prey cycles is the only source of instability as in the classical paradox of 
enrichment. The cycles pose a risk of extinction, occurin g stochastically in the low-point of the 



cycle or deterministically through subsequent bifurcations (|Andersenlll997l ). 



The onset of the predator-prey oscillations in the Hopf bifurcation relaxes the grazers control of 
the producers. Subsequently the density of producers grows, increasing the intraspecific competition 
between producers again while decreasing their nutrient content. At this point the effect of variable 
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conversion efficiency becomes important and tfius has to be taken into account in models. 

Due to tfie effect of intraspecific competition, further enrichment causes departure from the 
unstable region. The exact nature of the dynamics of the departure depends on the relation between 
food quality and competition. If intraspecific competition increases strongly while the food quality 
is still high, the oscillations disappear in another Hopf bifurcation. If however, the nutrient content 
decreases before the stabilizing effect of competition sets in then the system undergoes a saddle- 
node bifurcation in which a new stable state appears and a subsequent homoclinic bifurcation in 
which the oscillations end. 

Increasing the energy input further causes the extinction of the grazer population when the 
growth of the producer population decreases the producers nutrient content below the basic demand 
of the grazers. The extinction of the grazer can occur in a smooth transcritical bifurcation or in 
a saddle- node bifurcation giving rise to an Allee effect. In general the saddle-node bifurcation is 
encountered when the producer's growth is nutrient limited, so that a decreasing number of grazers 
leads to an increasing number of producers but docs not increase their nutrient content. 

If nutrient enrichment instead of energy enrichment is considered, the system will go through 
the same sequence of bifurcations in reverse order. However, in this case the sequence will end in 
the unstable region as further increase in nutrients will in general not impose constraints on energy, 
whereas in the energy-enrichment scenario increased energy supply induces a shortage of nutrients. 

The generic scenario discussed above shows that variable conversion efficiency has to be in- 
cluded in models to capture the dynamical constraints arising from stoichiometry when the grazer 
affects the biomass conversion efficiency. We emphasize that notable changes in the dynamics can 
already occur while the grazers nutrient limitation is still relatively mild. 

In the light of the generalized model the classical paradox of enrichment and the paradox of 
competition proposed here appear as two aspects of a more general paradox of constraints: Although 
constraints on the primary production are intuitively felt to be detrimental may benefit the system 
by lending stability. 

Apart from the specific results on stoichiometric consumer-resource systems the present work 
shows that it is advantageous to combine generalized and specific modeling. Generalized models 
avoid restricting functions in the model to specific functional forms. They can therefore be used 
to analyze certain dynamical properties of whole classes of specific models. In the present work 
we have shown that a single generalized model can accommodate the results of different specific 
models, which have been previously proposed. Thereby, generalized models can provide a unifying 
perspective, highlighting the differences and commonalities among different specific models. Con- 
versely, specific models provide more detailed insights and are necessary to study non-stationary 
dynamics that is not accessible in the generalized models. Because of this complementarity, we 
believe that the combination of insights from generalized and specific models will in the future 
prove to be a powerful strategy for the analysis of many complex ecological systems. 
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A. Online Appendix: Supporting material of the generalized modeling 

In this appendix we show some additional analysis of the generalized model. A more detailed 
description of the specific models can be found in the second appendix B. 



A.l. Normalization 

In the following we derive the normalized form of the generalized model corresponding to 

TtX = S{X,Y)-F{X)Y , 

= E{X,Y)F{X)Y - DY . ^ ' 

We assume that a positive steady state exist {X* > 0,y* > 0). For instance F* = F{X*) denotes 
the grazing rate in the steady state under consideration. As a first step of our analysis we define 
the normalized functions and variables shown in Tab. lAli By substituting these definitions into 
the original equations we obtain 

= ^i^s*s{x,y)-F*Y*f{x)y) , 
Ay = ^{E*F*Y*e{x,y)f{x)y-DY*y) . ^ ^ 

Let us now consider the system in the steady state {X* ,Y*) which corresponds to {x*,y*) = (1, 1) 
in the normalized variables. In the steady state the left hand side of Eq. ()A2p is zero and the 
normalized functions are by definition s(l, 1) = /(I) = e(l, 1) = 1. Therefore we obtain 

S* _ F*Y* 

E*F* _ D_ {^'^) 
Y* — y* ) 

which is reasonable as it states that the in the steady state the gain and loss rate of the pro- 
ducer (grazer) are identical. To simplify the equations we define := S* /X* = F*Y*/X* and 
ay := E*F* = D. The two quantities ax and ay are constants and can therefore be considered 
as parameters characterizing the steady state. From their definition it can be seen that these pa- 
rameters denote the biomass turnover rate of the consumer and grazer respectively. Using these 
parameters we can write the model as 



^x = ax{s{x,y) - f{x)y) , 
= oty{e{x,y)f{x)y -y) . 



(A4) 



As the final step we renormalize the timescale by a factor l/ox- In the normalized units of time 
the turnover rate of the producer is one whereas the turnover rate of the grazer is r := ay /ax- The 
new parameter r can therefore be interpreted as the biomass turnover of the grazer, measured in 
multiples of the biomass turnover of the producer. Writing the model in the rescaled units of time 
yields 



'^-^x = s{x,y) - f{x)y , 
^y = r{e{x,y)f{x)y -y), 



which are the normalized equations given in the paper. 
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Table 1. Bifurcation parameters of the generalized model 



Name Range Remarks 



r 




grazer turnover-rate 


< r < 1 


1 close to producer turnover-rate 


7 




sensitivity to producer 


> 


close to zero for saturated F{X), 

1 for F{X) linear in X 

2 for F{X) quadratic in X 


fx, 


ay 


sensitivity of production 
to producer and grazer 


ax < 1; cTy < 


substituted by and Cy 




Cy 


intra- and inter- 


< < 1 


no competition {S{X, Y) linear in X 






specific competition 


< Cj^ < 1 


and independent of y), 

1 only competition(S'(X*,y*) 0) 


Vx, 


Vy 


sens, of conv. efficiency 
to producer and grazer 


% < 0; % < 


substituted by and Uy 




Uy 


food quality 


< n:r < 1 
<ny <1 


1 for good food quality (constant E{X,Y)), 
for low food quality {E{X*,Y*) 0) 



Table Al. Normalized variables and functions 



Definition Substitution 



x:=^ X^X*x, 
y:=^ Y*y, 

s{x, y) := ^(Wly) s{X, Y) ^ S*s{x, y), 

e{x, y) := ^12^^ E{X, Y) ^ E*e{x, y), 

fix) := F{X) ^ F*e{x), 
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A. 2. Bifurcation conditions 

The Hopf bifurcation condition as well as the saddle-node bifurcation condition depend on the 
eigenvalues of the Jacobian ([5]). The saddle- node bifurcation occurs when a single real eigenvalue 
crosses the origin of the complex plane, acquiring a positive real part. In the moment of crossing 
the eigenvalue is zero and hence also the product of all eigenvalues, i.e., the determinant of the 
matrix must be zero. Hence, det(J) = is a necessary condition for the saddle-node bifurcation 
leading to Eq. 

In a Hopf bifurcation a complex conjugate pair of eigenvalues acquires positive real parts by 
crossing the imaginary axis of the complex plane. In the bifurcation the sum of the two eigenvalues 
vanishes. Because in the two-dimensional system, the sum of the two eigenvalues is the trace of 
the Jacobian, the condition trace(J) = is necessary for the Hopf bifurcation leading to Eq. ([9]). 
Furthermore, in the Hopf bifurcation the two eigenvalues crossing the axis are purely imaginary 
and of opposite sign. The product of these eigenvalues must therefore be a positive number. 
Thus, det(J) > is another necessary condition for the Hopf bifurcation leading to Eq. (llOp . We 
emphasize that for systems with more than two dynamical variables, computing the trace of the 
Jacobian is not sufficient for locating the Hopf bifurcation. However, even in this case explicit 
computation of the eigenvalues of the Jacobian is not necessary, as the Hopf bifurcation conditions 
may be directly obtained by the method described in (jGuckenheimer et al.lll997l : iGross and Feudel 



2004 ). 



A. 3. Supporting bifurcation diagrams 

In Sec. [3] we have stated that the specific value of the interspecific competition parameter Cy 
and the specific coupling between the food quality Ux and Uy do not change our results qualitatively. 
In this Section we support these statements by additional bifurcation diagrams. 

First, let us first discuss the effect of the inter-specific competition parameter Cy. As stated in 
Sec. [3] this parameter has no influence on the location of Hopf bifurcation points. Also, the effect of 
Cy on the saddle-node bifurcation surface in the variable-efficiency model is relatively minor unless 
Cy is comparable to Cx- As illustrated in Fig. Al (A, B) changing the value of Cy shifts the line in 
which the tangent bifurcation enters the parameter volume, but otherwise has little impact on the 
shape of the bifurcation surface. 

Second, we investigate the effect of changing and Uy independently. We compute a bifur- 
cation diagram for low competition (c^ = 0.01 and Cy = 0) and take and Uy as independent 
bifurcation parameters. As shown in Fig. Al (C), a decrease of either parameter, Ux or Uy, leads to 
an increase of the critical value of 7 at which the tangent bifurcation occurs. In addition increasing 
the parameter Ux also increases the values of 7 at which the Hopf bifurcation occurs. However, for 
moderate deviations from Ux = Uy the bifurcation surfaces remain qualitatively similar to Fig. 1. 
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Fig. Al. — Bifurcation diagrams of a generalized producer-grazer model. Hopf (red) and tangent 
bifurcation surfaces (blue) are shown. The fixed parameters are r = 0.3 (all), Cy = 0.6 (A), Cy = 0.95 
(B) and Cx = 0.01, Cy = (C). In all diagrams the steady state is only stable in the top volume. 
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A. 4. Replacement of the Hopf bifurcation as the primary source of instability 

In Sec. [3] we have shown visually that the Hopf bifurcation is replaced by a saddle-node bifur- 
cation when food quality is decreased. In the present section we support this visual impression by 
mathematical arguments. First, we show that no saddle-node bifurcations can occur if the efficiency 
is constant {ux = Uy = 1, rjx = rjy = Qi) and the functional response is monotonously increasing 
(7 > 0) before we proof that no Hopf bifurcation can occur if the food quality parameter is low, 
i.e., rix < 0.5 which is related to rjx < —1. 

The critical value of the grazer sensitivity where the saddle- node bifurcation occurs, 73, can 
derived from Eq. ^ as 

ay-l + r]y 

In models with constant conversion efficiency {r]x = rjy = 0) this equation simplifies to 73 = 0. 
Consequently, no saddle-node bifurcations can be observed in constant efficiency models of the 
form Eq. (lAip when the grazing rate increases monotonously with of the number of producers. 



From Eq. ([9]) and Eq. (|10p it follows that Hopf bifurcations can be observed at the critical 
value 7 = 7h where 

IH = -rrjy + (Jx, (A7) 

only if 

7H > 7S- (A8) 

This conditions shows that if a Hopf bifurcation exists, then it must be located at higher values of 
7 than the saddle-node bifurcation. 

Let us now assume that < 0.5 or equivalently r]x < —1. Since < 1 this implies 

{(Tx + Tlx) < 0. (A9) 
We multiply by the negative expression (ay — 1) and obtain 

{ay-l){ax + r]x)>0- (AlO) 
On the left hand side we add the positive term rrjyUy — rrjy + rrjy which yields 

r'qyOy - rr]y + rr]y + ay{ax + ijx) - {ax + r]x) > 0. (All) 
Finally we divide by the negative term ay — 1 + r]y and get 



rr]yay - rrjy + rr/^ + ayax + ayijx - ax - rjx 



< 0. (A12) 



ay-l + r]y 

Using Eq. (jATj) and Eq. ()A6P this is equivalent to 

7H - 7S < (A13) 
which contradicts the condition Eq. (jASp . Consequently, no Hopf bifurcations is possible if Ux < 0.5. 
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A. 5. The transcritical bifurcation 

The transcritical bifurcations in Sec. H] occur for c^; — )• 1. In the generahzed model, this 
situation is related to S{X* ,Y*) — > which is only possible if y* ^ as we can see from Eqs. ([T]- 
[2]) . Ecologically, this means that the producer becomes purely selflimited as the primary production 
and the grazer density Y* become zero. 

Mathematically, this case is problematic because positive, and therefore non-zero, densities 
have to be assumed in the derivation of the generalized model. Nevertheless, the bifurcation can 
be detected by considering the limit S{X*,Y*) — )• in which the normalization is strictly valid. In 
the limit ax = F*Y*/X* approaches zero, so that in the limit the condition for the saddle-node 
bifurcation 

det(J) = a^ayiriyiax - 7) - K " '^)iVx + 7))) = 0. (A14) 



is fulfilled. For a more rigorous discussion we refer the reader to IVan Voorn et al.l (|2008l ) 



A. 6. Higher codimension bifurcations and global dynamics 

At the intersections of bifurcation surfaces bifurcations of higher codimension are formed. 
These bifurcations can provide additional information about the system dynamics. The end of the 
Hopf bifurcation surface located in t he range 0.5 < < 1 is formed by a line of codimension-2 



Takens-Bogdanov bifurcation points (jKuznetsovl |2004| ) . While this bifurcation cannot be detected 



directly in experiments or observations it has important implications for the global dynamics. It 
is known that in addition to the Hopf and saddle-node bifurcation surfaces an additional surface 
of homoclinic bifurcation points emerges from the Takens-Bogdanov bifurcation line. Being a 
global bifurcation the homoclinic bifurcation cannot be detected directly in the generalized model. 
Nevertheless, the presence of the Takens-Bogdanov bifurcation already indicates that homoclinic 
bifurcations exist in models with variable conversion efficiency. 



B. Online Appendix: Supporting material of the specific modeling 



In this appendix we present the specific models from Sec. U] in greater detail, including the 
parameters that were used in our numerical investigations. For the sake of comparison we follow 
the notation of the original models, except when this would lead to unnecessary confusion. Note 
that small and capital letters are therefore no longer used to distinguish between normalized or 
non-normalized variables or processes. For examples of how the conven tional parar n eters used in 
specific models relate to the generalized parameters we refer the reader to lGross et al.l (120041 ) . where 
the case of the functional response 7 is discussed in detail. 
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B.l. DEB model and SDEB model 



In the DEB model the producer consists of two compartments. Assimilated nutrients are added 
first to a reserve or storage compartment, and then, in a second step, utilized for growth. Since the 
producers take up nutrients from the environment fast and efficiently, we assume that all nutrients 
are either in the structure or reserves of the producers X or in the structure of the grazers Y . For 
the the system formulation we follow iKooijman et al.l (|2007l ). neglecting the maintenance costs. 
Here the producer's reserve density m^v is obtained from the conservation of nutrient in the system 



mN{X, Y) = N/X - TINY Y/X - n^x 



(Bl) 



for a total constant amount of nutrient in the system. The chemical indices njsfx and n^y denote 
the producers' and the grazers' nutrient content per carbon. This implies that X{t) G {0, N/njyx) 
and P{t) G [0,N/nNp). 



Following iMuller et al.l (120011 ) . the processes in Eqs. ()Aip are given for the DEB model by 

kNmNiX,Y) 



S{X,Y) 



VNX + mN{X,Y) 



nX) = f^ and 
E{X,Y) = (y-^ + {yYNmN{X,Y))-^ 



{VYX + VYNmNiX, Y)y 



(B2) 

(B3) 
(B4) 



where the growth rate of the producers S{X, Y) is assumed to follow Droop-kinetics and the specific 
feeding rate F(X) is the Holling-type-II functional response. The conversion efficiency E{X, Y) 
of the grazers resul t s from the SU r ules for the parallel processing of complementary compounds 
(jO'Neill et al .11 19891 : lKooiimanll2000l ). The parameters yyx and yyN respectively denote the yield 
of the producer's structure and reserves when converted to grazer growth. 

For the simplified constant efficiency version of the DEB model we assume that the grazer 
consumes the structural part of the producer only. Then the growth rate of the grazer follows a 
Holling-type-II functional response, such that 



F{X) 
E{X,Y) 



K + X' 

VYX- 



(B5) 
(B6) 



The parameters of both models are summarized in Table IBll 



LKE model 



The LKE model was proposed and analyzed in (iLoladze et al.l l2000l ) . The model explicitly 
focuses on the essential nutrients carbon and phosphorus. In contrast to the DEB models no 
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storage of nutrients is represented directly in the model. The density of phosphorus ri(X, Y) in 
the producers population X is variable but not less than a minimal density q. The density of 
phosphorus 9 in the grazer population Y is assumed to be constant. The primary production 
follows a logistic growth. If the system is limited by energy input then the carrying capacity is 
assumed to be a constant K. If however phosphorus is limiting then the carrying capacity is given 
by the upper limit for the producer density, i.e., the total available phosphorus {P — 9Y) divided by 
the minimal phosphorus density in the primary producer q. Hence, the classical carrying capacity 
in the logistic growth is replaced by the minimum function min(i^, (P — 9Y)/q). 

The grazer consumes the producer's carbon at the rate F{X), where F{X) is the functional 
response. At the same time the producer's phosphorus is consumed at the rate r]{X, Y)F{Y). If the 
growth of the grazer is carbon-limited the conversion efficiency is assumed to be a constant E. But, 
if the phosphorus density of the producer r]{X, Y) is below the phosphorus density of the grazer 
then the conversion efficiency is decreased by the ratio r]{X,Y)/9. Consequently the conversion 
efficiency is defined as E(X,Y) = mm{E,Er]{X,Y)/9). In summary, the processes in Eqs. ()A1|1 
are in the LKE model given by 

- P - 9Y 

E{X,Y) = mm{E,E7]{X,Y)/9) with r]{X,Y) = — (B7) 

X 

a + X 

SiX,Y) = bx(l- ■ , ) (B9) 

V mm{K/X,r]{X,Y)/q) J 

where the constant P denotes the total amount of phosphorus in the closed system. Note that the 
grazer egests nutrients that are not used for growth. The egested products are mineralized and 
sequestered by the producer instantaneously. As a result, no external carbon and phosphorus pools 
are assumed. 

In the model the efficiency E(X,Y) is constant if carbon is limiting, but becomes inversely 
proportional to X if phosphorus is limiting (cf. Eq. ()B7P ). In the generalized parameter space 
this is related to a switch of Hx from 1 to 0.5, i.e., r]x switches from to 1. As we have shown in 
Sec. IA.4l the value Hx = 0.5 is exactly the threshold where the Hopf bifurcation vanishes. In order 
to yisualize the switch a two-parameter bifurcation diagram of the specific model, similar to Fig. 5 



in (ILoladze et al.ll2000l ). is shown in Fig. Bl (C). In the figure the parameter values at which the 
switch takes place is marked by the curve T. Additionally a homoclinic bifurcation curve, G, two 
transcritical bifurcation curves, TCi, TC2, a Hopf bifurcation H, and a saddle-node bifurcation 
S, are shown in the figure. We note that the Hopf bifurcation as well as the homoclinic and the 
saddle-node bifurcation end at the switch. In contrast to the results from the generalized model 
the Hopf and saddle-node bifurcations appear on the same side of the switch in the specific model. 
To understand why this difference arises note that the bifurcations take place on different steady 
states, and hence correspond to different values of the generalized parameters. For the steady state 
undergoing the Hopf bifurcation the switch T is from = 1 to n^, = 0.5. By contrast, for the 
steady states undergoing the saddle- node bifurcation the switch T is from n^; = 0.5 to n^j, = 1. 
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Tablc Bl. Parameter table of the DEB model 





Name 


Value 


Units 


N 


Total nutrient in the system 


0-8.0 


mol 1"^ 


riNX 


Chemical index of nutrient in X 


0.25 


mol mol-i 


uny 


Chemical index of nutrient in Y 


0.15 


mol mol-i 


Un 


Reserve turnover-rate 


0.25 


h-i 


Vnx 


Yield of N on X 


0.15 


mol mol-i 


jPm 


Maximum specific assimilation rate 


0.4 


mol mo\~^ 


K 


Half saturation constant 


10 


mM 


Vyx 


Yield of Y on X 


0.5 


mol mol-i 


Vyn 


Yield of Y on N 


0.8 


mol mol-i 


jxAm 


Maximum specific assimilation rate 


0.15 


mol mol-i 


D 


Hazard rate of grazer 


0.005 


-1 



Table B2. Parameter table of the LKE model 





Nam(> 




Uuitri 


p 


Total phosphorus 


0.025 


mg P 1-1 


E 


Maximal production efficiency in carbon terms 


0.8 




b 


Maximal growth rate of the producer 


1.2 


day-i 


D 


Grazer loss rate (includes respiration) 


0.25-0.27 


day-i 


e 


Grazer constant P/C 


0.03 


(mg P)/(mg C) 


q 


Producer minimal P/C 


0.0038 


(mg P)/(mg C) 


c 


Maximum ingestion rate of the grazer 


0.81 


day-i 


a 


Half-saturation of grazer ingestion response 


0.25 


mg C 1-1 


K 


Producer carrying capacity limited by light 


0.25-2.0 


mg C 1-1 
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In the generalized model we observe that a decrease of competition always tends to destabilize 
the steady state. At low food quality where no Hopf bifurcations can be found this destabilization 
is caused by a saddle-node bifurcation. We see from Fig. Bl (C) that, in the LKE model the same 
bifurcation scenario can be found when phosphorus is limiting: an increasing total phosphorus 
concentration can lead to a destabilization due to the saddle-node bifurcation S. 

While the presence of the homoclinic bifurcation confirms our expectations from the generalized 
model, the Takens-Bogdanov bifurcation from which the homoclinic bifurcation emerges in the 
generalized model, does not appear in Fig. Bl (C). This difference arises because the switch of 
rix from 1 to 0.5 avoids the parameter region where the Takens-Bogdanov bifurcation is located. 
Numerical continuation of the homoclinic connection terminates at the switch T. Therefore, one 
could ask whether the homoclinic bifurcation can still be linked to the Takens-Bogdanov bifurcation 
of the generalized model. To confirm that the Takens-Bogdanov bifurcation is indeed the organizing 
center of the Hopf and homoclinic bifurcations a smooth analogon to the LKE model, which avoids 
the discontinuity, is studied in the subsequent section. 



Smooth analogon model 



To overcome discontinuities we formulate a new model (SA), constituting a smooth approxi- 
mation of the LKE model. The producer is assumed to consists of two components: a phosphorus 
pool and the structure, which has a fixed stoichiometry given by the P/C ratio q. 



To model t he assimilation of the producer (Eq. ()B10p the SU- formulation (jO'Neill et al 



198S 



Kooijmanll20ld ) is used where both nutrients, carbon and phosphorus are assumed to be essential. 
In analogy to the LKE model we assume absence of phosphorus in the environment. The total 
phosphorus density in the producer is the same as in the LKE model, i.e., rj^X, Y) = (P — 9Y)/X. 
Hence the phosphorus of the pool is r]{X, Y) — q. In the LKE model the light energy is represented 
by a light-limited carrying capacity K. By contrast we assume an external carbon pool C instead 
representing the energy resource for the producer. Consequently, the growth of the producers 
depends on carbon influx from the environment proportional to C — x and internal phosphorus from 
the pool P — 9y — qx. Let us emphasize that the assumption for the carbon influx is necessary to get 
a good agreement with the LKE model. We note that this model is not closed since carbon is further 
converted into grazer biomass. However, the limitation of the carbon flux can be interpreted as a 
simple formulation of limited photosynthe tic capacity due to self -shading of producers. The system 
is closed for nutrients but open for energy (jKooiiman et al.ll2002l ). For obtaining an analogon to the 
LKE model, we assume that light is not limittin g. Then, the carr ying capacity can be interpreted 
as the total amount of carbon, C, in the system (jKooi et al.lll998l ). 
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r 
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Fig. Bl Onc-paramctcr and two-paramctcr bifurcation diagrams of tlic LKE mode 1 iLoladze et al.ll2000l ) (A, C) and the smootli analogon 

model (B, D). The upper panels (A,B) represent bifurcation diagrams corresponding to the line P — 0.025 in the lower panels (C.D). In both models 
along the line P — 0.025 a positive, stable steady state emerges from a transcritical bifurcation TCi , (solid line) and becomes unstable (dashed line) 
at the Hopf bifurcation H. A stable limit cycle that emerges from the Hopf bifurcation vanishes in a saddle-node homoclinic bifurcation G. Here a 
pair of steady states, one stable and one unstable emerge from the saddle-node bifurcation, S. The stable solution of the saddle-node bifurcation 
exchanges stability with the zero equilibrium in another transcritical bifurcation TC2 ■ The lower diagrams show that the saddle-node homoclinic 
bifurcation G turns in both models into a saddle homoclinic bifurcation for lower values of P. In (C) the curve T marks the stoichiometric switch 
of the grazer minimum function, i.e. ((P — By) / x) / — 1. The Hopf bifurcation curve H , the homoclinic bifurcation curve, G, and the saddle-node 
bifurcation curve, S, all terminate in the curve T. By contrast, the diagram for the smooth analogon (D) shows that both the Hopf bifurcation 
curve H and the homoclinic bifurcation curve G terminate in a Takens-Bogdanov bifurcation point TB. The saddle- node bifurcation Si terminates 
together with another saddle-node bifurcation S2 in a cusp bifurcation point N. 
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The processes of the SA model in Eqs. ()Aip are given by 

S{X,Y)=bX J. J. ^-^ (BIO) 

^ ^ 1 _|_ ^PC J ^PC A PC ^ ' 

^ {C-X)Bc (P-eY-qX)Bp {C-X)Bc + {P-eY-qX)Bp 

. T , Kpc , Kpc Kpc 
with := 1 + — — + 



CBc PBp CBc + PBp 

E(X,Y)=E r^TTTT- (Bll) 

^ + TKj^x) ~ ^'^^ s — ) 

so that b is the maximum initial producer-growth-rate. The parameter 9 is again the phosphorus 
density in the grazer and r]{X, Y) the phosphorus density in the producer Eq. (jB7|l . The parameters 
Be and Bp denote the assimilation preferences of the producer for C and P respectively. The 
parameter Kpc is a saturation constant. 

The consumed amount of carbon and phosphorus by the grazer are proportional to F{X) while 
ri(X, Y)F{X) respectively. Note that there is no distinction between phosphorus originating from 
the structure of the producer's structure and phosphorus pool. While the use of the SU-formulation 
in the derivation of Eq. (jBlip requires the fluxes to be independent, the application of this formalism 
is justified by assuming that after ingestion both nutrients from the assimilation (catabolic) process 
become available for growth as unrelated chemical substances while both being essential. We have 
chosen the parameters of the SA model (see Table IB3|) such that a good correspondence to the 
LKE model is achieved. 

The bifurcation diagram of the smooth SU-model formulation shown in Fig. Bl (B) is very 
similar to the results from the LKE model (Fig. Bl (A)), where the total carbon concentration C 
now takes the place of the carrying capacity K. Again the appearance of the saddle-node and the 
homoclinic bifurcation are in agreement with the results from the generalized analysis. 

From the generalized analysis, we expect a Takens-Bogdanov bifurcation to be the organiz- 
ing center of the Hopf and the homoclinic bifurcations. However, the saddle-node and the Hopf 
bifurcations cannot meet in a Takens-Bogdanov bifurcation since both belong to different steady 

Table B3. Parameter table of the smooth analogon to the LKE model 





Name 






Value 


Units 


E 


Yield of carbon and phosphorus 






0.96 




Kpc 


Saturation constant 






1 


mg c 


Be 


Producer assimilation preferences 


for 


C 


0.002 


1 (mg C)-i 


Bp 


Producer assimilation preferences 


for 


P 


2 


1 (mg P)-i 



Note. — 6, C, P, 0, ry, c and e are the same as in table [B2] 



- 33 - 



states. We therefore need to find a saddle- node bifurcation of the steady state undergoing the Hopf 
bifurcation. Indeed by increasing D shghtly from 0.25 to 0.27 we observe that the steady state that 
becomes unstable in the Hopf bifurcation undergoes a saddle-node bifurcation S2 and turns into a 
stable steady state again in the saddle- node bifurcation 5*1. The resulting bifurcation diagram is 
shown in Fig. 3 (top right). 

From the generalized analysis we expect that the Hopf bifurcation H and the saddle-node 
bifurcation point 5*2 to meet in a Takens-Bogdanov bifurcation, from which also the homoclinic 
bifurcation emerges. A two-parameter continuation of the bifurcations shown in Fig. Bl (D) con- 
firms this expectation. Decreasing the total phosphorus content P causes the Hopf bifurcation line 
to terminate in a Takens-Bogdanov bifurcation. Prom this point also the homoclinic bifurcation in 
which the limit cycle is destroyed emerges (see Fig. Bl (B)). Furthermore the diagram shows that 
the saddle-node bifurcation S2 emerges together with Si from a cusp bifurcation N. At d = 0.25 
we are above the ^2 curve and therefore the second bifurcation is absent in Fig. Bl (B). 

In summary, compared to LKE model the SU formulation yields qualitatively similar results. 
However, the continuous model allows us to link the avoidance of the paradox of enrichment to 
an underlying Takens-Bogdanov point which was concealed by the discontinuous behavior of the 
original model. 



Model by Diehl 2007 



A model proposed by Diehj ( 2007 ) considers a producer population X and a grazers population 



y in a mixed water column of depth z. It is assumed that the producer assimilates the available 
nutrient R fast and efficiently until the algal nutrient content (quota), per carbon Q, reaches a 
certain maximum Qmax- The specific algal growth rate p is assumed to depend on the quota Q and 
on the local light intensity. The latter is described in dependence of the local depth s and X by the 
Lambert-Beer's law I{X, s) = li^e"^^^'^^^''^^^ where k is the producer light attenuation coefficient 
and Khg the background light attenuation coefficient. It is assumed that inorganic carbon is never 
limiting and the contribution of the grazer on shading is neglected. Respiration is included for both 
species, with a rate Im for the producer and a rate m for the grazer. The processes of the model 
can be written as 

z 

S{X, y) = § J PiHX, s),Q{X, Y)) ds - ImX, (B12) 


F{X) = (B13) 



O(X.Y) ^1 



Q(X, Y) = miniQmax, ?i2L_SL), (b15) 
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The functional response F{X) is of Holling-type-II. Similarly to the smooth analogon model, the 
form of the variable conversion efficiency, E[X,Y), is derived by the synthesiz ing units ap proach 
and is equivalent to the expression in Eq. (IBTT]) and the expression given in (jPiehll 120071 ). Here 
the growth of the grazer is limited by algal carbon and the nutrient quota Q{X, Y) of the con- 
sumed producer biomass. The assimilation rate of carbon is assumed to be 50% (c=0.5) while the 
assimilation rate of nutrients is assumed to be 100%. 

The local specific algal growth rate p{I{X, s)^Q{X^Y)) is assumed to be co-limited by the 
algal nutrient-to-carbon ratio Q{X^Y) and the local light intensity I{X,s). This co-limitation is 
modeled by a product of two monotonously increasing, saturating functions of Q(X, Y) and I{X, s) 
respectively, 

I{X,S) Qr. 



p{I{X,s),Q{X,Y))=p^^^^^^ ^^^^ 
which leads by integration to the primary production rate 



1 



QiX,Y) 



S{X,Y) 



1 Pn 



zkA + Kb 



In 



H + z) 



1 



Qr 



Q{X,Y) 



ImX. 



(B16) 



(B17) 



Finally, note that the DLL model above was studied by iDiehll (j2007l ) as a limit case of a 4- 
dimensional system which models dissolved and sedimented nutrients explicitly. Differences between 
both model arise mainly for shallow water columns due to the sedimentation of algal nutrients. 
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Table B4. Parameter table of the DLL model 







Value 


Units 
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D 
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day-^ 
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LO 


day-i 


H 


Half-saturation constant for light-dependent 
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120 mmol photons m~^ s~^ 
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